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I. INTRODUCTION 



A. THE IDEA OF ARMA MODELING 

The goal of linear modeling is to accurately represent an observed data sequence 
as the output of a linear filter. The idea of representing a complicated process with a 
comparatively simpler model has many different applications. Curve fitting in math- 
ematical modeling, analysis of electronic devices using equivalent circuits, and system 
transfer functions in automatic control are just a few examples outside the area of dig- 
ital signal processing. Parametric modeling has also a large number of applications in 
signal processing. Currently there is considerable interest in the parametric modeling 
approach to spectral estimation. In speech processing the applications include digital 
transmission, storage, and synthesis of the speech signal. Our particular interest is 
the modehng of sonar signals, such as biologies and other underwater acoustic data. 
This work forms part of an overall research program in sonar signal modeling. The 
research will help to understand the relative benefits of signal domain algorithms 
versus algorithms based on coefficients of the transfer function. It is hoped that the 
method developed in this thesis will become an important tool in the overall effort 
for sonar signal modeling. 

In linear modeling the filter used to generate the data sequence is usually rep- 
resented by a linear difference equation with constant coefficients. The 2:-transform 
of this type of system is a rational polynomial function. Three type of models are 
derived from this kind of systems; they are known as autoregresive (AR), moving 
average (MA), and autoregresive moving average (ARMA). 

Much work has been done on AR models, which correspond to all-pole systems. 
The reason for that is that the parameters for the model can be obtained by solving 
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linear equations, and a great body of theory has been developed that applies to 
this problem [Ref. 1, 2], Relatively much less work has been done with MA and 
■ARMA models. However, since MA models have limited applications and are almost 
as difficult to obtain as ARMA models, most interest centers on the latter. The 
filter in this type of models has both poles and zeros, and the design fundamentally 
involves nonlinear equations. A properly designed ARMA model can provide better 
performance than an AR model, with a smaller number of parameters. ARMA 
modeling is the topic of this thesis. 

B. WHY AN ITERATIVE ALGORITHM IN THE SIG- 
NAL DOMAIN 

There are many different approaches to the problem of ARMA modeling. The 
majority of them are based on statistical techniques [Ref. 3j. Some of these methods 
regard the data as a realization of a random process while others focus on the data as 
given [Ref. 1]. Data oriented methods try to minimize some criterion that estimates 
how well the model fits the data, in most cases the least squares error between the 
signal and the model. Stochastic approaches may attempt to estimate the model 
parameters directly from the data by solving nonlinear equations or by spectral fac- 
torization. The maximum likelihood procedure, for example [Ref. 1, 4], is essentially 
nonlinear. A number of indirect methods have been developed that modify the norm 
of the error by separating the AR and MA parts of the problem so that at least some 
of the equations to estimate the parameters become linear. This approach is found 
in procedures such as the Prony’s method. Shank’s method, and the least squares 
modified Yule- Walker method [Ref. 1, 2, 5, 6]. A different type of approach replaces 
the nonlinear problem with iteration while trying to solve for the AR and MA pa- 
rameters simultaneously. The iterative prefiltering method of Steiglitz and McBride 
[Ref. 7, 8] is of this type. 
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Our method is also of the latter type. However, the advantage is that it works 
directly with the poles of the rational model, which we know affect the performance 
of the system. The poles are displaced in specific directions so that the new model 
minimizes the error between the model output and the original signal. Iterative 
prefiltering on the other hand works with the coefficients of the transfer function, 
so it is difficult or impossible to predict its effects on the poles and zeros of the 
system. The new algorithm is much more dependable with respect to convergence 
than the iterative prefiltering algorithm because it moves poles and zeros specifically 
to minimize the error between the model and the original signal. 

C. THESIS OUTLINE 

The remainder of this thesis is organized as follows. Chapter II presents the 
modeling methods that are used in this thesis and gives a brief explanation of all 
of them. Chapter III introduces the reader to the theory of multidimensional op- 
timization by gradient methods and develops the iterative Prony method, which is 
the main contribution of this thesis. Chapter IV presents the results of testing the 
algorithm on simulated and real acoustic data and compares these results with those 
obtained using iterative prefiltering. Chapter V gives conclusions and suggestions for 
future research. 
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II. ARMA MODELING OF SIGNALS 



A. MODELING METHODS USED IN THIS THESIS 



This thesis deals with deterministic approaches to ARAIA modeling. Two types 
of modeling methods are considered. First we have non-iterative methods hke Prony 's 
method and its alternate signal domain form [Ref. 1]; second we consider iterative 
methods like iterative prefiltering and the new iterative Prony method developed in 
this thesis. 

The goal of Prony ’s method is to represent a given sequence x[n] as the impulse 
response of a linear time invariant (LTI) system. In the transform {z) domain this 
representation has the form 



X{z) 



B{z) 

A(z) 



( 2 . 1 ) 



where X{z) is the z-transform of x[n] and B{z)/ A{z) represents the transfer function 
of the system. This approximation is explained in more detail in the next section. 
What has become known as Prony’s method in the current signal processing literature 
differs from Prony’s original work in some respects. Our basic form of Prony’s method 
solves for the coefficients of the transfer function in (2.1). 

The signal domain form of Prony’s method is closer to Prony’s original work 
[Ref. 1, p.560] and seeks to represent the data in terms of a set of damped exponen- 
tials as 



j;[n] Cir^ + C 2 T 2 H h Cprp 



where the r^ are the roots of the denominator polynomial A{z) and the Ck are the 
complex coefficients required for the expansion. Both forms involve linear equations 
and least squares techniques. 
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~l ;0 e[n] = x[n] — x[n\ 



x[n\ 



B(-A Ztohz-'' 

■4(2) 1 + E£=1<H2”‘ 



x[n 



Figure 2.1: Block diagram for the direct method for signal modeling. 

The iterative prefiltering method used is due to Steiglitz and McBride [Ref. 7, 8] 
and attempts to match the data with a model of finite order using an efficient iterative 
approach. It differs from Prony’s method in that it solves for both numerator and 
denominator polynomial coefficients simultaneously at each iteration. 

Whenever a signal is modeled using a fixed-order rational polynomial model, an 
approximation has to be made and some kind of measure has to be used to determine 
the “goodness” of the model. In this thesis the least squares error norm [mboxl 2 
norm) is used to measure the approximation error. This norm measures the energy 
of the error and is the norm most widely used primarily due to its mathematical 
tractability [Ref. 2]. 

B. OVERVIEW OF MODELING METHODS 
1. Prony’s Method 

The derivations of all the modeling methods presented in this section follow 
those in [Ref. 1, pp. 550-564] and [Ref. 2]. As stated above, Prony’s method aims 
at representing a signal as the impulse response of an LTI system. Figure 2.1 shows 
an implementation of this approximation which is known as the direct method. The 
system function X(z) in the transform domain is a rational polynomial function 
B{z)fA(z) with Q zeros and P poles. The error signal e\n\ is computed as the 
difference between the response of the system x[n] and the given signal x[n], i.e. 

e\n] = x[n] — x[n]. (2.2) 
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The LTI system is chosen to minimize the sum of squared errors 



<^D = ^|eN|^ (2.3) 

n 

This problem leads to nonlinear equations whose solution (if a unique solution actu- 
ally exists) turns out to be a very difficult task [Ref. 1.2]. To avoid these difficulties, 
a number of indirect methods have been developed for modeling. Prony’s method is 
one of these procedures and can be derived as follows. The LTI system satisfies the 
difference equation 

.tfn] -f aix[n — 1] -(-••• + api’[n — P] = 6o6[n] -f- 6i6[n — 1] -[-••• + — Q\. 

(2.4) 



If the requirement that 



ir[n] = x[n], n = 0, 1, . . . , iV^ — 1, 



is applied to (2.4), where Ns is the length of the data, and the difference equation is 
evaluated for n = 0, 1, . . . , A^s — 1, the result is the matrix equation 



■ x[ 0 ] 


0 




0 


... 0 






■ 6 o • 


x[l] 


x[ 0 ] 




0 


... 0 






bi 


x[2] 


x[l] 




x[ 0 ] 


... 0 




■ 1 




62 


• 






• 


... ; 




«i 






: 








... 1 




02 




j 


x[Q\ 


x[Q - 


1] 


x[Q - 2 ] 


••• x[Q-P] 






— 


bq 


x[Q 1] 


x[(5] 




x[Q - 1] 


— P + 1] 




dp 




0 


. x[Ns - 1] 


x[iV,- 


- 2 ] 


x[Ns - 3] 


• • • x[Ns — P — 1] ] 






. 0 



This can be written as 



Xb 

X.4 



a = 



b 

0 



( 2 . 6 ) 
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where 



Xb = 



■ x[0] 


0 


0 


... 0 


x[l] 


x-[0] 


0 


... 0 


x[2] 


x[l] 


x[0] 


... 0 


. ^[Q] 


x[Q - 1] 


1 

... 14 


X[Q-P]_ 



and 



X. = 



x[Q + 1] x[Q] 



x[Q - 1] 



(2.7) 



x[Q — P + 1] 



( 2 . 8 ) 



x[Ns — 1] a:[A^s — 2] x[Ns — 3] • • • x[Ns — P — 1] 
and a and b are the vectors of coefficients appearing in (2.5). The lower partition 



X^a = 0 (2.9) 

represents an overdetermined set of linear equations that need not have an exact 
solution. This set of equations can be solved by least squares, where a is chosen to 
minimize the least squares norm of the equation error = ||e^||^ in 



X^a = e^. 



( 2 . 10 ) 



This leads to a set of linear equations called the normal equations, which can be 
written compactly as [Ref. 1, pp. 536-537] 




( 2 . 11 ) 



and can be solved for a and Sa- Once the vector a is known, it can be substituted 
in the upper partition of (2.6) 

b = Xga (2.12) 



to solve for the vector b. Although it is referred here simply as Prony’s method, 
this procedure is also known as the modern Prony method or the extended Prony 
method. 
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e^[n] = x[n] * a 



i,...,6q,0.0.0.. 



[n] - b[n] 



Figure 2.2: Block diagram for the indirect modeling problem. 

Although Prony’s method is simple to implement, it is important to keep 
in mind that it is an indirect method. In particular this procedure (see Fig. 2.2) 
minimizes the squared magnitude of the error 



e, 4 [n] = a:[n] * a[n] — 6[n] 



(2.13) 



where the sequences a[n] and b[n] are defined as 



def 



and 



r 1 

a[n\ = 



6[n] 



a„; 0 < n < P (gq = 1) 

0; otherwise 



a„; Q <n<Q (gq = 1) 
0; otherwise 



(2.14) 



(2.15) 



where P and Q are the number of poles and zeros respectively. Equation 2.13 repre- 
sents a different least squares error from that in (2.3) where the quantity to minimize 
was the squared magnitude of the error e[n] (see Fig. 2.1). The practical significance 
of this difference is that frequently there is a loss of accuracy in estimating the poles 
and zeros by Prony's method [Ref. 3]. 



2. Signal Domain Form of Prony’s Method 



An alternative formulation of the method described above can be obtained 
if the problem (2.1) is stated in the signal domain by representing a:[n] in terms of a 
set of complex exponentials: 

x\n\ ~ Cjr" C 2 V 2 + • ■ • cpv^ (2.16) 
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where as stated before, the are the roots of the polynomial A{z), assumed to be 
distinct, and the complex coefficients c* provide for the linear combination of the P 
roots. 

This approximation can be initially formulated as in the section above and 
(2.9) can be solved in the least squares sense for the coefficients of A{z) (ie. the 
vector a). The roots r^. of A{z) can then be found and (2.16) can be evaluated for 
n = 0, 1, . . . , — 1 to produce the set of equations 



“ 1 1 


•• 1 








■ a;[0] 


r\ V2 


• • rp 

. . r2 
Tp 




C2 


= 


^[1] 

^[2] 


1^1 ^2 


f— 1 
1 

CQ 

. . 




. CP . 




x[Ns — 1] 



This set of equations can then be solved in a least squares sense to obtain the vector 
of coefficients c. 

In the case of multiple roots at the same location, a slight variation of the 
same procedure can be used. Suppose, for example, that r\ is a double root. In this 
case, the approximation is 



a;[n] ^ Cir" + C 2 nr" + h cpVp 



(2.18) 



and the matrix equation to solve for the coefficients becomes 



1 



0 

2r\ 



1 

rz 



^N,-l 

'1 



. N ,-\ N ,-\ 



1 

rp 



^N,-\ 









1 

1 








x[l] 




^2 


= 


x[2] 




. cp . 




x[Ns — 1] 



(2.19) 



(JV,-l)rr ,3 ... 

This situation is rare, however, because computational errors and errors inherent to 
the modeling method itself contribute to produce roots that may be very close to 
each other, but not exactly at the same location. 
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3. Iterative Prefiltering 



The iterative prefiltering method attempts to solve the '‘direct” problem 
mentioned in subsection 1 of this chapter and to refine the initial pole-zero estimate 
by solving a succession of linear problems. Equation 2.2 (error for the direct problem) 
can be written in the ^-domain as 



E[z) = X{z) - X{z) = X{z) - 



B{z) X{z)A{z) - B(z) 



( 2 . 20 ) 



A(z) A{z) 

The notion of iteration can be introduced that allows computation of a new set of 
poles and zeros based on the last known set of poles. Iterative prefiltering replaces 
the error for the direct problem at the (i + 1)‘^ iteration with the iterative error 
function 

If is the inverse 2-transform of \jA^''\z), then the error can be written in the 

signal domain as 



= a:[n] * (2.22) 



The coefficients and 6b+i)|j2] are selected to minimize the corresponding sum 

of squared errors 

5(.+i) ^ ^ |e<‘+^)[n]| (2.23) 

n=P 

at each iteration; this situation is shown in Fig. 2.3. No general proof of convergence 
has been given for this algorithm; however, it is easy to see that if the iteration does 
converge, it must produce the same answer as the direct method. Specifically, at 
convergence and (2.21) becomes the same as (2.20). 

If we use an indirect modeling procedure like Prony’s method to compute 
the initial vector a and we define as 



a:[n] * /i**^[n]. 



(2.24) 
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Figure 2.3: Block diagram for the iterative prefiltering method. 



then the sequences and for n = 0, 1, . . . , Ns — 1 can be computed from 

the recursive difference equations 

P 

/ib)[n] = ^ — k] (2.25) 

k=i 



and 

p 

^ (2.26) 
fc=i 

Thus the error (2.22) can be written as 

e(’+i)[n] = ^ - j]. (2.27) 

k=0 j=0 



In order to find the coefficients and the error is written in matrix form 

for n = 0, 1, . . . , — 1 as 



Xb) 



ab+i) ■ 

_bb+i) 



- eb'+i) 

'I 



where 



= 


a;(d(pj 

x^^P+1] 


x^^^[P-l] ■■■ 

xb)[P] 


a;(*)[0] 

xb)[l] 




. xW[A^s - 1] 


xb)[7V5_2] ... 


.rb)[A^5 - 1- P] _ 


= 


1] 


fi(*^[P-l] ... 


h(i)[P-Q^l] 




. - 1] 


/i(*)[fV5-2] ... 


h(i)[Ns-l-Q] . 



(2.28) 



(2.29) 



(2.30) 
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and 





■ 1 




[ 1 


= 










. '4*" . 




. r" . 



(2.31) 



Equation 2.28 is analogous to (2.10). Thus the least squares problem defined by 
minimizing the norm of the error in (2.28) can be reduced to 



X(0 H(’) ] Hh) 



a(i+i) 



5h+i) 

0 

0 



(2.32) 



where 



5(1+1) ^ ||efi+‘>f . (2.33) 

This linear equation can then be solved for the v^ectors of filter parameters. 



At this point, it is instructive to compare the performance of the three 
modeling methods outlined in this chapter by applying each to model a small segment 
of a transient sound corresponding to a wrench being struck. This wrench sound was 
recorded and sampled in the laboratory at a sampling rate of approximately 10,240 
Hz. This signal, denoted by wrenOl, was also used and modeled in [Ref. 9]. Figure 
2.4 shows a 100 point segment of the signal wrenOl with one of the three models 
(Prony’s, signal domain form of Prony’s, and Iterative Prefiltering ) overlaid. In all 
three cases the signal was modeled with four poles and four zeros. As can be seen, 
the difference between iterative prefiltering and the first two models is significant. 
The non-iterative methods match only the initial points of the sequence, and produce 
poor approximations in modeling the remaining part of the signal. This is due, in 
large part, to poles not sufficiently close to the unit circle. Iterative prefiltering, 
on the other hand, produces a model which is close to the real sequence along the 
complete segment. 
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magnitude magnitude magnitude 



2000 






Figure 2.4: Signal wrenOl and its 4 poles/4 zeros models, (a) Using Prony’s method, 
(b) Using Signal Domain form of Prony’s method, (c) Using Iterative Prefiltering. 
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III. ITERATIVE ALGORITHM IN THE 
SIGNAL DOMAIN 



A. MULTIDIMENSIONAL OPTIMIZATION BY GRADI- 
ENT METHODS 

A multidimensional function f{x\,X 2 .,---,Xn) that is continuous and differen- 
tiable can be minimized using one of several very powerful hillclimbing techniques 
known as gradient methods [Ref. 10, p. 84]. Some of those methods are derived on 
the basis of a quadratic model that can be obtained from a truncated Taylor series 
expansion of /(x). Let denote the value of x at the k'’^ iteration. Then for any 
point X = x^*^ -f 6 ; when S is small, the function can be approximated by 

-f 5) « -f (3.1) 

where g and G represent the vector of first derivatives and the matrix of second 
derivatives of the function /(x) respectively and they should be available at every 
point. In Newton's method the iterate is taken to be x^*^ -f where the 

correction 6^'^^ minimizes This method is only well defined when the matrix 

of second derivatives G is positive definite, in which case the k'’^ iteration of Newton’s 
method is given by the following procedure [Ref. 11, pp 44-46]: 

1. solve = — g^*’^ for 6 =6^^^ 

2. set = x^^’^ -f 

(3.2) 

The fact that G^^^ may not be positive definite when x^^'^ is far from the solution, 
and that even when G^^l is positive definite convergence may not occur, makes this 
method undesirable as a general formulation of a minimization algorithm. However, a 
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number of variations to the basic method have been proposed that are more suitable 
for a general class of problems. One of these methods is Newton's method with line 
search [Ref. 11, pp. 47-49] in which the Newton algorithm is used to generate a 
direction of search 

(3.3) 

which can later be used in a line search algorithm to actually calculate the correction 
S. In the cases when is not positive definite, the linear search can be made 

along choosing the correct sign to ensure a descent direction. However, some 

difficulties that arise here (like very high numerical costs and failure of convergence 
for some special cases) make this an undesirable approach for our algorithm. 

As stated before, Newton's method is defined only when the matrix is 

positive definite, and this matrix is positive definite only when the error 6 is “small”; 
or better stated, the method is defined only in some neighborhood of in 

which ^^*^^(<5) agrees with 4- ^) in some sense. In such cases, it is correct to 

choose + with the correction minimizing for all x^^l + 3 

in This method is referred to as the restricted step method because the step is 

restricted by the region of validity of the Taylor series. [Ref. 11] 

The region of definition for the iteration can be expressed as 

= {x ; [jx - x^'^)]] < (3.4) 

where jj • [j denotes the norm of the vector. In this case, the optimization problem 
can be stated as; 

minimize^ subject to [|<5[[ < (3-5) 

As mentioned before, the least squares norm is the one most commonly used in this 
type of problems, so it is the one used in this thesis and is denoted as || • [[ 2 . The 
problem that now becomes apparent is how to select the error margin of the 
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neighborhood (3.4). This margin should be as large as possible because the iteration 
step is directly related to it. Various methods have been proposed to control the 
parameter one of these methods attempts to insure that the Newton’s search 

direction problem (3.3) is always defined [Ref. 11, pp. 100-103]. It does so by adding 
a multiple of the unit matrix I to and computing the new problem 

(g<‘I + ^l) = -g<‘l (3.6) 

where the net effect is that increases in u cause ||<^||2 to decrease, and vice versa. 

If we define 

/(/:) 

then the ratio represents a measure of the accuracy to which approxi- 
mates -h on the step, and as the accuracy increases gets closer to 

unity. Using (3.7), Marquardt [Ref. 12] suggests an algorithm that tries to adaptively 
maintain as large as possible while controlling the ratio The iteration 
of such an algorithm is stated as: 

1. given and calculate and 

2. factor -f if not positive definite, reset and repeat; 

3. solve (3.6) to find 

4. evaluate /(x^^^ and hence 

5. if < 0.25 set 

else if r^^'> > 0.75 set /2 

else set j/U+0 = 

6. if < 0 set = x^*^ else set x^*'*'*' = x^*) 4- 

Here the parameters 0.25, 0.75, 4 and 2 are arbitrary, and > 0 is also chosen 
arbitrarily [Ref. 11, pp. 102-103]. Proofs of global and second order convergence for 
this algorithm are given in [Ref. 11. pp. 96-98] for the cases when the first and sec- 
ond derivatives of the function /(x) exist, and the vector x^*^ belongs to a bounded 
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7 i-dimensional space for all k. Although this method does have some disadvantages, 
it represents a good basis for the formulation of a general minimization algorithm. 
Some variations of this method were considered, but it was found that for the spe- 
cific application of .ARMA modeling in the time domain, these variations produced 
an extremely high overhead in calculations. 

B. THE ITERATIVE PRONY METHOD 

Let us now return to the problem of representing a sequence 2 r[n] as a linear 
combination of complex exponentials. Equation 2.17 can be written as 

Rc = X -I- € (.3.8) 



where e is called the equation error, x represents the data which may or may not 
be complex, c is the vector of complex coefficients, and R is the matrix of complex 
roots, which can be written more specifically as 







1 






1 






1 








^Ri 


+ 


jrii 


rR2 


+ 


;>/2 


TRp 


+ 


JT'Ip 




R = 


i^Ri 


+ 


jrR? 




+ 


jri2? 


■ • • (’’flp 


+ 


jri,? 


(3.9) 




_ 


+ 




(^H2 


+ 




••• i'^Rp 


+ 


jrip)^‘-^ 





where and r/_ represent the real and imaginary components of the root, re- 
spectively, and Ng is the number of data samples. 

By defining 

Q ||€||2 = €*^e = (Rc — x)*^ (Rc — x) , (3.10) 

it is clear that the problem is to find the vector r = i^R + jri of P complex roots that 
minimizes the function Q(r). 

The first and second derivatives of Q with respect to r are represented by the 
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vector g and matrix G, which are defined by 



and 



G = 



r >90 

do 



r)0 



dvR 

dtf 

dri 

dd 



dQ 



drj 



P -I 



(3.11) 



d'^Q 


a^Q 


d^Q 


a'^Q 


d^Q 


a^Q 1 


a^a 


dVR^ dVR^ 

d^Q 


drRj^ drRp 

d^Q 


drR dri^ 
d^Q 


drR drj^ 
d^Q 


drR drip 

d^Q 




drR^drR^ 


drR^drRp 


arR^drj^ 


arR^arj^ 


avR^arjp 


a'^Q 


a^Q 


a^Q 


a^Q 


a^Q 


920 


dr R parity 

a^Q 


drRpdrR^ 

d^Q 


drRpdrRp 

d^Q 


arRpari^ 

a^Q 


drRpdrj^ 

d^Q 


arRpdrjp 
92 0 


dri^drR^ 

a^a 


drR^ 

a^Q 


drj^ drRp 
d^Q 


drj^ drj^ 
d^Q 


drj'drj^ 

d^Q 


dri^arip 

a^Q 


drj^ arR^ 


drj^dvR^ 


drj^dTRp 


dri^arj^ 


ari 2 arj^ 


drj^arjp 


a^Q 


d^Q 


a^Q 


a^e 


92 Q 


a^Q 


^ dripdvR^ 


drip dr R^ 


drjpdrRp 


drjpdri^ 


drjpdrj^ 


dripdrip ^ 



(3.12) 



In order to provide for a more compact notation, define the gradient operator with 

respect to the complex vector r as the vector of partial derivatives 

a 



Vr = 



V 

Vr 



rR 



dr 






drR^ 



dVRp 

d 

dri^ 

d 

drr 



drj 



(3.13) 
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consequently, (3.11) can be expressed as 

g = VrQ = 

Equation 3.12 can also be written as 






(3.14) 



G = 



a 


\ do 


dQ 


dQ 


dQ 


dQ 


dQ 


a 


' dQ 


drR^ 

dQ 


dVRp 

dQ 


dri^ 

dQ 


dri^ 

dQ 


drip 

dQ 




dvR^ 


drR^ 


dVRp 


drj^ 


dri^ 


drip 


a 


\ dQ 


dQ 


dQ 


dQ 


dQ 


dQ 


drRp 




drR.^ 


dVRp 


ar/j 


drj^ 


drip 


a 


' dQ 


dQ . . 


dQ 


dQ 


dQ 


dQ 


ar/j 

a 


drR^ 
' dQ 


dVR^ 

dQ 


dTRp 

dQ 


drj^ 

do 


drj^ 

dQ 


drip 

dQ 


drj^ 




drR^ 


dTRp 


drj^ 


dri^ 


drip 


a 


dQ 


dQ 


dQ 


dQ 


do 


dQ 


dT,p 


dVR.^ 


dVR^ 


dTRp 


drj^ 


drj^ 


drip 



(3.15) 



Now using a somewhat more convenient notation, (3.15) can be rewritten as 



G = 



d 


■ dQ 


dQ ■ 




dVR^ 


drR^ 


dri^ 




d 




dQ ' 






drR^ 


dri^ 





i = \ ... P 
k=\...P 



(3.16) 



and it becomes clear that the matrix of second derivatives can be expressed compactly 



as 



G = Vr (VrQ) 



= Vr 



(Vr„Q)’^ (V„Qf 



(3.17) 



The following subsections derive explicit expressions for the two quantities g and G. 
1. Vector g of first derivatives 
From (3.8) it is easy to see that 

(3.18) 



de 

dr; 



OR 



-t:—c and — — 
dri dvi 



= c 



dri 
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where r,- represents the real or imaginary parts of the root. Using (3.18) and the 
chain rule in (3.10) leads to 



dQ 5R .7^ 

= e*-' c + C ^ 



5R* 



-e , 



drn, drn, 

and explicit evaluation of the partial derivatives of the matrix R results in 



(3.19) 



dQ 

drn. 



= 6*^ 



0 

1 

2{rn, +;>/,) 
^{rR,+jrif 



Cl + 



cj 



. (iV, - l)(rH, ^ 

0 1 2 (cr, - jr/. ) 3 (cr, - jrj, f 



(3.20) 



where Cj is the component of the vector c. If the quantity is defined as 



0 

1 



def 

^ t 



2(rR. +jrj,) 

3(rR, +jriy 

(A'. 



(3.21) 



then 



dR 




dvR,"" 






= 


dvR, 



(3.22) 



(3.23) 



and (3.20) can be written compactly as 



,t*T 

= e ^t+^, €• 



dr 



R, 



(3.24) 
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At this point it can be recognized that (3.24) represents the addition of two scalar 
quantities where the first one is the complex conjugate of the second, so (3.24) can 
be further simplified to 



do 

drn. 



= 2 Re 



ii « 



(3.25) 



where Re[-] denotes the real part of the vector. Finally, using (3.25) for f = 1 • • • P 
in the upper partition of (3.14) produces 



= 2 Re < 





r 1 


X 




^ 1 






S2 




< 




e ^ 




/>v 


> 



(3.26) 



A similar procedure can be used to obtain the gradient of Q with respect 
to the imaginary part of the vector r. Once more, from (3.18) and the chain rule 
apphed to (3.10), it follows that 



drj^ 

drj^ 



c = 



= -jCi- 



(3.27) 



(3.28) 



These results can be used to generate an expression similar to that in (3.24) for the 
vector of first derivatives of Q with respect to the imaginary components of r 

dQ 



dr 7 






(3.29) 



which in turn defines the vector of partial derivatives with respect to the imaginary 
components as 



Vf; Q = 2 Im < 






L^: 



T -I 



. 



(3.30) 
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Equations 3.26 and 3.30 can now be combined as shown in (3.14) to obtain the final 
expression for the vector of first derivatives 



g 



^v,Q 



= 0 



Re < 



T -I 



Im < 






L 4p J 



■■ 



• P -I 



(3.31) 



2. Matrix G of second derivatives 

An expression for the matrix G of second derivatives can be obtained as 
follows. Substituting (3.24) and (3.29) into (3.16) yields 



G = 





j 


.•"s.-sfe] ) ■ 


. 4: ( A, 


j 





i = l...P 
k=l...P. 



(3.32) 



Then using the chain rule and both expressions in (3.18) leads to 
G = 









arR, arR, 



dri. ^ dr I « 1 ^ « I dr, ^ dTj. ^ 



i = \ . . . P\ k = I ... P. 



•T _ a£^ 

drji ' drrt. dro. Qtc 






• r 

dvR ^ dTR ^ 



drr^ 



. *T 



^ drr^^t drr^^ 



(3.33) 



From the definition of ^ in (3.21) it is seen that 






dr 



i c / 

- = Sik = 

Rk 



0 

0 

2 

6(rn. +;V/.) 
12 (rn, +jriy 



_ (yV,-l)(iV,-2)(rH.+jr;.) 



Ns-3 



0, 



In the same way it can be shown that 

dif 

9rR, 



= S 



ik 



c^^ if i=k 



otherwise. 



(3.34) 



(3.35) 



drr 



= JSik 



(3.36) 



drn 



■ *T 
-J^ik ■ 



(3.37) 



These last four expressions together with (3.22), (3. 23), (3. 27), and (3.28) can now be 
substituted in (3.33) to obtain 



G = 









*T^ I t*Te ^*T 

e Sik + ik ~ ^k ~ SjJt e 






i = I ... P 
k = i...p. 



(3.38) 



Finally notice again that the elements of the matrix G are formed by addi- 
tions and subtractions of complex scalars with their respective complex conjugates; 
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thus the final expression for G is given by 



2 Re 




+ 2 Re [s*Je 


2 Im 


if it 


-2 Im 






—2 Im 




+ 2 Im [s'Je 


2 Re 




- 2 Re 


[srJ.; 





i=l...P: i = (3.39) 

where it should be remembered that 8,^. = 0 for all i ^ fc. 

3. Algorithm implementation 

An iterative method for ARMA modeling in the time domain was imple- 
mented using the results of the last two subsections in conjunction with the algo- 
rithm presented in section A of this chapter. We call this method the iterative Prony 
method. 

The algorithm uses the signal domain form of Prony’s method to calculate 
an initial model for the given sequence. From there it uses the calculated model 
to compute the error e, the vector of first derivatives g, and the matrix of second 
derivatives G and iterates until specific conditions are met. Figure 3.1 is an example 
of how the algorithm changes the position of the poles and zeros of the initial model 
in order to minimize the error. This figure represents the poles and zeros of a transfer 
function of order (4,3) — 4 poles 3 zeros — that was overmodeled using a (6,5) order 
model. It is clear from the figure that the tendency in this case is to have a pole-zero 
cancellation (see second and third quadrants) of two poles and two zeros as expected. 

Some features were added to the basic algorithm in order to deal with special 
modeling cases. Specifically, if the initial model has some roots on the real axis, then 
because of the way the algorithm iterates, those roots never move away from the real 
axis. A modification was therefore introduced to deliberately displace those roots 
from the real axis and proceed with the iterations. If the tendency of the roots is 
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1.5 



1 - 
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I 
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Figure 3.1: Displacement of the poles and zeros of an iterative Prony’s model 

“to go back” to the real axis, then they are returned to their initial position, and the 
iterations continue. Otherwise the roots may continue to spread apart and move as 
a complex pair. Figure 3.2 is an example of the displacement of the poles and zeros 
of an order (4,3) model. In this case the modeled signal actually has two poles on 
the real axis, and the initial model correctly placed two of the poles in the real axis. 
Those poles are displaced from the real axis by the algorithm, but then after some 
iterations it is clear that the poles are tending to return to the real axis. At this point 
the poles are forced back to the real axis by setting their imaginary parts to zero 
and the iterations continue until convergence is obtained. The opposite situation is 
shown in Figure 3.3. In this case the initial model also has two poles located on the 
real axis, but contrary to the case presented above, the roots, after being displaced 
from the real axis, continue to move away from the axis until they reach their final 
position in the first and fourth quadrants closer to the unit circle. 
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Figure 3.2: Displacement of the poles and zeros of a 4-3 model. The modeled signal 
in this case actually has two poles on the real axis 

1 

0.8 
. 0.6 
0.4 
0.2 

I ^ 

- 0.2 

- 0.4 

- 0.6 

- 0.8 

-1 



* = poles before iterations 




- 


X = poles during iterations 


X 

X 


f 


0 = zeros during iterations 


X 

-i 








_ 








> o o* 0 o o o oooooo 


* 


- 




' - ^ 




* 


n 


f 






X 

X 




- 




» 


- 




Re[] 



* = poles before iterations 


- 


X = poles during iterations 


o 


0 = zeros during iterations 


s 


A 


V' 

* 

•DO HO i - 
\ ^ 

T 




X 

Red 



Figure 3.3: Displacement of the poles and zeros of a 4-3 model. The initial model 
shows two poles in the real axis, the final model in this case has no poles in the axis 



IV. PERFORMANCE OF THE ALGORITHM 
AND MODELING RESULTS 



A. TEST DATA USED IN THIS THESIS 

Two types of signals were used in this thesis to test the performance of the 
algorithm. The first type, which will be called simulated test data, consists of five 
sequences (tOl to t05) each one hundred points long that were produced as the 
impulse response of a known rational system. Noise to produce an SNR in the range 
of 10 to 15 dB was added to the original sequences, and the resulting sequences were 
designated as t01_n to t05_n. The original signals are described in Table 4.1 by their 
transfer functions and the location of their poles and zeros. 

The second group of test signals consists of recorded acoustic data. Two of these 
signals were recorded and sampled in the laboratory. One of them is the sequence 
wrenOl already mentioned in Chapter II; the other one wais obtained from human 
speech, in particular, the signal voweLa corresponds to 100 samples of the Spanish 
vowel a. The remaining three signals from the group of acoustic data were recorded 
at sea by a submarine platform; they correspond to sounds produced by marine life 
and ice cracking. The description of the acoustic signals is presented in Table 4.2. 

B. PRESENTATION OF RESULTS 

All simulated test signals were modeled twice with the iterative Prony method. 
In the first test the exact number of poles and zeros of the original model was used; 
in the second test all signals were modeled using two more poles and zeros than the 
original model. This last test is considered closer to a real life situation where the 
exact order of the signal to be modeled is unknown. 
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TABLE 4.1; DESCRIPTION OF SIMULATED TEST DATA 



NAME 


TRANSFER FUNCTION 


POLES 


ZEROS 


tOl 


IJ(^\ 1-0.772“^ 

— 1-1.5392-^+0.9052-2 


0.9513 Z± 0.6286 


0; 0.770 


t02 


IJ(^\ — 0.7442-^-1.69932-2+0.6302“^ 


0.9422 L ± 0.6200 
0.8987 Z± 0.6385 


0; 0.4686 
1.8069; oo 


^ 1-2.9772-1 +3.9092-2 -2.5202-'i+0.7172-^ 


t03 


_ 0.17U-1-0.155Z-3 


0.9512 Z ± 1.8846 
0.9514 Z ± 2.3041 


0; 0.9521 
—0.9521; oo 


V / 1+1.8612-1+2.5882-2 + 1.6842-^+0.8192-1 


t04 


fJ(^\ — 0.98l2“^-1.2792“2+0.87l2“^ 


0.9513 Z ± 0.2097 
0.9049 Z ± 2.0943 


0: oo 

0.9423 Z ± 0.8068 


l-0.9562-i+0.0402-2-0.7052-'i+0.74l2-i 


t05 


U ( ^\ 1—0.752“^ 

— 1-1.782-1+0.792-2 


0.8441; 0.9358 


0; 0.750 



TABLE 4.2: DESCRIPTION OF ACOUSTICTEST DATA 



SIGNAL 


DESCRIPTION 


wrenOl 


Transient sound corresponding to a wrench being struck 


voweLa 


Spanish vowel a 


bio_2133a 


Sperm whale 


bio_2385a 


Porpoise whistle 


bio_80a 


Ice cracking 
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To mathematically produce a meaningful measure of the performance of a mod- 
eling algorithm can be quite difficult since various norms can be deceiving when 
comparing errors of signals with large differences in magnitude. Two different ap- 
proaches to measure the performance of the algorithms were therefore used in this 
thesis. The first is quantitative and involves computing the squared-norm of the error 
between the model and the actual signal and dividing it by the total energy (norm) 
of the signal. The second approach is to overlay in a plot the model and the original 
signal in order to provide a visual comparison of the results. This is less quantitative 
but frequently more revealing of errors in the modeling process. 

1. Simulated test data 

The first data sets modeled were the simulated test data sets. Figure 4.1(a) 
is a comparison of the normalized errors that result when the sequence t01_n is 
modeled with 2 poles and 1 zero using both iterative prefiltering and iterative Prony 
methods. Figure 4.1(b) and (c) show 100 points of the sequence t01_n and the 
two order (2,1) models. At this point there is no noticeable difference between the 
iterative prefiltering and the iterative Prony models. Figure 4.2(a) again shows a 
comparison of the normalized errors between an iterative prefiltering model and an 
iterative Prony model of t01_n for the case when the signal t01_n was overmodeled 
using models of order (4,3). Although the difference between the two modeled signals 
in this case is not large, notice that the error for the iterative prefiltering method 
initially increases before decreasing while the error for the iterative Prony method 
decreases monotonically. This is the first example of a pattern that repeats in all but 
one of the simulated test signals that were modeled. Figures 4.3 through 4.10 give 
similar comparisons for the remaining simulated signals. The pattern, which can be 
seen in Figures 4.1 to 4.10, is that when the signals are modeled with a number of poles 
and zeros different from that of the actual order of the signal (always overmodeling 
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Figure 4.1: Signal tOl-ii and its 2 poles-1 zero models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal tOl-n and 
the iterative prefiltering model, (c) Signal tOl-u and the iterative Prony model. 
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Figure 4.2: Signal tOl-ti and its 4 poles-3 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal tOl-n and 
the iterative prefiltering model, (c) Signal tOl-U and the iterative Prony model. 



31 



magnitude magnitude Error 






Figure 4.3: Signal t02-n and its 4 poles-3 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t02~n and 
the iterative prefiltering model, (c) Signal t02.n and the iterative Prony model. 
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Figure 4.4: Signal t02-n and its 6 poles-5 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t02-n and 
the iterative prefiltering model, (c) Signal t02-n and the iterative Prony model. 
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Figure 4.5; Signal t03.n and its 4 poles-3 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t03-n and 
the iterative prefiltering model, (c) Signal t03-n and the iterative Prony model. 
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Figure 4.6: Signal tOSjn and its 6 poles-5 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t03^n and 
the iterative prefiltering model, (c) Signal t03-n and the iterative Prony model. 
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Figure 4.7: Signal t 04 ~n and its 4 poles-3 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal 104 -n and 
the iterative prefiltering model, (c) Signal t04~n and the iterative Prony model. 
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Figure 4.8: Signal t04~n and its 6 poles-5 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t04~n and 
the iterative prefiltering model, (c) Signal t04-n and the iterative Prony model. 
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Figure 4.9: Signal t05-n and its 2 poles-1 zero models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t05-n and 
the iterative prefiltering model, (c) Signal t05-n and the iterative Prony model. 
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Figure 4.10: Signal t05-n and its 4 poles-3 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal t05-n and 
the iterative prefiltering model, (c) Signal t05-n and the iterative Prony model. 
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in this case), the behavior of the iterative prefiltering method tends to degrade (i.e., 
it takes longer for the algorithm to reach convergence) while the behavior of the 
iterative Prony method remains the same or even improves as in the case of t05_n 
shown in Figures 4.9(a) and 4.10(a). Parts (b) and (c) of Figures 4.1 through 4.10 
show the original signals and their respective models overlaid. It can be seen that 
the models arrived at by both methods follow the original signals very closely in all 
cases. Table 4.3 lists the location of the poles and zeros of the systems used to model 
all the simulated noisy signals. These systems were obtained using both the iterative 
prefiltering and the iterative Prony methods. The poles and zeros shown in Table 4.3 
can be compared to the poles and zeros of the original simulated signals presented in 
Table 4.1. It is clear that the location of the poles and zeros of the modeled signals 
should not be exactly the same as those of the original signals because some noise (in 
the order of 10 to 15 dB SNR) was intentionally added before the modeling process. 
However, Tables 4.1 and 4.3 show a close relation between the location of the poles 
and zeros of the original sequences and the position of the poles and zeros of the 
modeled signals. 

2. Acoustic test data 

As mentioned in the last section the acoustic test data represents sounds 
recorded both underwater and in a laboratory environment. In some cases shorter 
segments were selected for modeling due to the complexity of these signals. Once 
again the iterative prefiltering algorithm was used, and its results were compared to 
the results obtained using the iterative Prony method. 

Before presentation of results, it is important to explain how the model 
produced by the iterative prefiltering method was selected. For all the cases, the 
same number of iterations was used both for the iterative Prony and for the iterative 
prefiltering methods. In order to obtain the best possible results from the iterative 
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TABLE 4.3: POLE-ZERO LOCATION OF THE SYSTEMS USED TO MODEL 
THE SIMULATED NOISY TEST DATA 



SIGNAL 
NAME k 
(ORDER) 


ITERATIVE 

PRONY 


ITERATIVE 

PREFILTERING 


POLES 


ZEROS 


POLES 


ZEROS 


t01_n 

(2.1) 


0.9511 Z ±0.6291 


0; 0.7612 


0.9507 Z ± 0.6290 


0; 0.7619 


tOlji 

(4,3) 


0.9510 Z ±0.6290 
0.9217 Z± 3.1397 


0.7636 

0.9476: 0.9151 


0.9506 Z ± 0.6289 
0.9867 Z ± 2.0944 


0.7630 

0.9886 Z ± 2.0880 


t02ji 

(4,3) 


0.9545 Z ± 0.6238 
0.8455 Z ± 0.6743 


13.2102 
2.1568; 0.1799 


0.9250 Z± 0.6313 
0.9142 Z± 0.6238 


6.8483 

1.6619; 0.4650 


t02ji 

(6,5) 


0.9541 Z ± 0.6228 
0.8505 Z ± 0.6797 
0.9896 Z ± 2.7699 


17.1527 
2.1237; 0.0952 
0.8760 Z ± 2.7876 


0.9301 Z ± 0.6335 
0.9068 Z± 0.6173 
0.9221 Z ± 2.4366 


6.5260 

1.8203; 0.5542 
0.9800 Z ± 2.3896 


t03_n 

(4,3) 


0.9510 Z± 1.8850 
0.9505 Z ± 2.3029 


22.9950 
1.0511; 0.9365 


0.9499 Z± 1.8851 
0.9514 Z± 2.3030 


12.2468 
0.9680; 0.8020 


t03_n 

(6,5) 


0.9509 Z ± 1.8851 
0.9506 Z ± 2.3029 
0.9356 Z ± 0.2779 


21.1093 
0.9340; 0.6882 
1.1445 Z± 0.3507 


0.9492 Z ± 1.8854 
0.9515 Z± 2.3028 
0.9796 Z± 1.5986 


7.5747 

0.9566; 0.6099 
0.9870 Z± 1.6035 


t04_n 

(4,3) 


0.9022 Z ± 2.0981 
0.9515 Z± 0.2097 


22.9853 

0.9494 Z ±0.8212 


0.9048 Z ± 2.0956 
0.9516 Z± 0.2090 


2362.7 

0.9000 Z ± 0.8090 


t04_n 

(6,5) 


0.9040 Z ± 2.0966 
0.9517 Z ± 0.2094 
0.7237; 0.7221 


17.5349 

0.9404 Z ± 0.8333 
0.7099 Z± 3.1147 


0.9053 Z ±2.0933 
0.9518 Z± 0.2093 
0.9542; 0.8105 


57.6799 

0.9486 Z ± 0.8189 
0.9461; 0.7803 


t05_n 

(2.1) 


0.9416; 0.7840 


0.6693 


0.9363; 0.8325 


0.7236 


tOSji 

(4,3) 


0.9366; 0.8345 
0.9571 Z ± 2.8523 


0.7330 

0.9571 Z ± 2.8436 


0.9368; 0.8231 
0.8839 Z ± 2.7995 


0.7019 

0.8998 Z ± 2.7977 
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prefiltering algorithm (especially in those cases when convergence was not obtained), 
the model selected was that corresponding to the iteration with the lowest error 
between the model and the original sequence. Figure 4.11 is an example of one of 
the cases when convergence was not reached using the iterative prefiltering method. 
It can be seen that if the system selected is the one obtained at the last iteration 
(20‘^) of the algorithm (which in this case corresponds to a peak of the error), then 
the resulting model does not follows the original signal at all. However, if we select 
the system from the iteration for which the error is the smallest (17‘^ iteration in this 
case), then we obtain a model that is closer to the original signal. Some insight can 
also be obtained if we look at the behavior of the poles and zeros of the system while 
the error of the iterative prefiltering algorithm is oscillating. Figures 4.12 and 4.13 
show the position of the poles and zeros during the oscillation that exists between 
iterations 15 to 20 of the case presented in Figure 4.11. For the first part of these 
iterations when the error is low, the poles remain at approximately the same position. 
At iteration 20, (Fig. 4.12(f)) the poles suddenly spread apart, some to even outside 
the unit circle, causing the large increase in the error. A similar effect occurs with 
the zeros, although there is some significant movement of the zeros in the earher 
iterations. This type of behavior is avoided in the iterative Prony method where a 
controlled and systematic displacement of the poles and zeros is produced to finally 
reduce the error between the modeled and original signals. 

The same approach used in the last subsection was used to model the acous- 
tic data. All signals were first modeled with a low order system and subsequently 
remodeled using a higher order system. The results are shown in Figures 4.14 to 4.23. 
In most of the cases the models obtained using the iterative Prony method closely 
follow the original signals. It can also be observed that the iterative Prony method 
does not have as large a dependency in the order of the model selected as does the 
iterative prefiltering method. While it is of course necessary to select a model with a 
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Figure 4.11: Signal voweLa being modeled by iterative prefiltering using a (10,9) 
order system, (a) Normalized squared-norm of the error between the model and the 
actual signal, (b) Signal voweLa and the iterative prefiltering model selected from 
the 20‘* iteration, (c) Signal voweLa and the iterative prefiltering model selected 
from the 17‘^ iteration. 
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Figure 4.12: Behavior of the poles of the system used to model the signal voweLa 
during one oscillation of the iterative prefiltering algorithm, (a) iteration, (b) 
16'^ iteration, (c) iteration, (d) 18'^ iteration, (e) 19'^ iteration, (f) 20‘^ 
iteration. 
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Figure 4.13: Behavior of the zeros of the system used to model the signal voweLa 
during one oscillation of the iterative prefiltering algorithm, (a) 15*^ iteration, (b) 
16*^ iteration, (c) iteration, (d) 18*^ iteration, (e) iteration, (f) 
iteration. 
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high enough order to properly model the original signal, once such a model has been 
selected, increments or small variations to its order do not produce degradations on 
the performance of the algorithm. On the contrary, it was found that in some cases 
the performance of the iterative prefiltering algorithm can be significantly reduced 
when the order of the model is increased slightly. It can also be seen that the itera- 
tive Prony algorithm tends to provide a closer match to the data than the iterative 
prefiltering method, and in most of the cases the rate of conv'ergence of the iterative 
Prony method was higher than that of iterative prefiltering. Another important point 
that can be extracted from the results presented in Figures 4.14 to 4.23 is that while 
convergence with neither algorithm is guaranteed, we obtained conv^ergence with the 
iterative Prony method in aU cases for these acoustic signals. The same was not 
true for iterative prefiltering. In most of the cases the error for the iterative Prony 
method begins to decrease starting at the first iteration, and although the change is 
not monotonic in all cases, the error after a few iterations is consistently lower than 
the initial error. 
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Figure 4.14: Signal wrenOl and its 7 poles-6 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal wrenOl and 
the iterative prefiltering model, (c) Signal wrenOl and the iterative Prony model. 
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Figure 4.15: Signal wrenOl and its 12 poles-11 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
wrenOl and the iterative prefiltering model, (c) Signal wrenOl and the iterative 
Prony model. 
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Figure 4.16: Signal voweLa and its 10 poles-9 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
voweLa and the iterative prefiltering model, (c) Signal voweLa and the iterative 
Prony model. 
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Figure 4.17: Signal voweLa and its 14 poles-13 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
voweLa and the iterative prefiltering model, (c) Signal voweLa and the iterative 
Prony model. 
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Figure 4.18: Signal bio-2133a and its 8 poles-7 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
bio^2l33a and the iterative prefiltering model, (c) Signal bio-2133a and the iterative 
Prony model. 
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Figure 4.19; Signal bio-2133a and its 12 poles-11 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
bio-2133a and the iterative prefiltering model, (c) Signal bio-2133a and the iterative 
Prony model. 
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Figure 4.20: Signal bio.2385a and its 8 poles-7 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
bio^2385a and the iterative prefiltering model, (c) Signal bio-2385a and the iterative 
Prony model. 
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Figure 4.21: Signal bio-2385a and its 12 poles-11 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
bio-2385a and the iterative prefiltering model, (c) Signal bio-2385a and the iterative 
Prony model. 
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Figure 4.22: Signal bio. 80 a and its 8 poles-7 zeros models, (a) Normalized squared- 
norm of the error between the models and the actual signal, (b) Signal bioSOa and 
the iterative prefiltering model, (c) Signal bioSOa and the iterative Prony model. 
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Figure 4.23: Signal bioSOa and its 12 poles-11 zeros models, (a) Normalized 
squared-norm of the error between the models and the actual signal, (b) Signal 
bio-SOa and the iterative prefiltering model, (c) Signal bioSOa and the iterative 
Prony model. 
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V. CONCLUSIONS 



A. DISCUSSION OF RESULTS 

In this thesis, a new method for modeling signals in the time domain is developed 
and applied to model both recorded acoustic data and simulated signals produced 
as the impulse response of a known system. We call this method the iterative Prony 
method. In most of the simulated test data sets the models provided by the iterative 
Prony method are sufficiently close to the original signals; in most cases, it is diffi- 
cult to distinguish between the signal and the model. When modeling the acoustic 
data distortion becomes apparent in some of the models, which may be due to the 
complexity of the structure of the signals. However, this distortion is no worse than 
for any of the best algorithms that have been used to model this data previously. 

The new algorithm was compared very specifically to the iterative prefiltering 
algorithm [Ref. 7, 8] which has been used in modeling a variety of acoustic data 
[Ref. 3, 9]. The rate of convergence of the iterative Prony method was in most of the 
cases comparable or superior to that of the iterative prefiltering algorithm. Thus, 
while iterative prefiltering sometimes has convergence problems, the new algorithm 
is much more dependable in that respect. The price to pay for this improvement is 
in the number of computations. While the number of floating point operations per 
iteration in the iterative prefiltering method is approximately 

64(P + Q - 1)3 + 8A,(P H- Qf -h 10(P -f Q)N, -f 12N,, 
iterative Prony requires about 

672P3 -h (24fV, -b 102)p2 + (eOA^, -f- 46)P -b 198iV, 
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floating point operations at each iteration. For example for a complex data set of 
length 100 (iVj = 100) and P = Q — 8 we have approximately 452.400 floating 
point operations per iteration using iterative prefiltering versus 572,360 using iter- 
ative Prony. If we increase the order of the model to P = Q = 16, then we have 
approximately 2,787,824 operations per iteration in the iterative prefiltering algo- 
rithm versus 3,509,560 in the iterative Prony method. 

B. RECOMMENDATIONS FOR FUTURE WORK 

The iterative prefiltering algorithm has been the main tool in the modeling 
efforts for sonar data modeling [Ref. 9, 13]. The new iterative Prony algorithm is 
now at a stage where it can be substituted for the iterative prefiltering algorithm and 
tested in operational use. To do so needs some further programming to make the 
segmentation of the data automatic and to make the entire modeling procedure more 
of a “turn crank” operation. These should be some of the very next steps. In addition, 
the practical implications of the increased computation needs to be addressed, and 
if possible new methods need to be developed to help reduce computations. 

In a larger sense the work reported in this thesis can be used as a base for 
possible applications of the iterative Prony method in the problems of filter design, 
speech processing, and spectral estimation. The expressions for the vector of first 
derivatives g and the matrix of second derivatives G of the error derived in Chapter 
III can be used along with different minimization methods to provide for other new 
modeling methods that may adapt better to specific modeling problems. 
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